A guide for Biologists and Ecologists
15th July 2024
Associate Professor in Data Science and Genetics at the University of East Anglia.
Academic background in Behavioural Ecology, Genetics, and Insect Pest Control.
Teach Genetics, Programming, and Statistics
Norwich - “A Fine City!”
I hope you end up with more questions than answers after this workshop!
Reproducible
Robust
Transparent
Reusable
Shareable research materials
| dataset | mean_x | std_dev_x |
|---|---|---|
| away | 54.26610 | 16.76982 |
| bullseye | 54.26873 | 16.76924 |
| circle | 54.26732 | 16.76001 |
| dino | 54.26327 | 16.76514 |
| dots | 54.26030 | 16.76774 |
| h_lines | 54.26144 | 16.76590 |
| high_lines | 54.26881 | 16.76670 |
| slant_down | 54.26785 | 16.76676 |
| slant_up | 54.26588 | 16.76885 |
| star | 54.26734 | 16.76896 |
| v_lines | 54.26993 | 16.76996 |
| wide_lines | 54.26692 | 16.77000 |
| x_shape | 54.26015 | 16.76996 |
Originally created by Alberto Cairo in Download the Datasaurus: Never trust summary statistics alone; always visualize your data
Weissgerber TL, Milic NM, Winham SJ, Garovic VD (2015) Beyond Bar and Line Graphs: Time for a New Data Presentation Paradigm. PLoS Biol 13(4): e1002128.
Easy enough to rapidly prototype graphics at the “speed of thought”
Powerful enough for final “publication”
Also see ggbeeswarm
A model is a representation of what we think is going on in the real world. It can be a good fit, or a poor fit.
A well fitted model could be used to make predictions about what might happen in the real world.
A poorly fitted model could also make predictions but they might be disastrously wrong.
If we want our inferences to be accurate the model must represent the data collected.
{fig-alt=’ Model bridges - Andy Field Discovering Statistics’ width=485}
If data is symmetrically distributed, the mean and median will be close, especially as n increases.
This is also know as a normal distribution
Understanding the shape of our data informs the summary statistics we can use
Depending on our question, if the data is not symmetrical this could really matter…
Box bounded by 25th and 75th percentile (first and third quartile). This is called the interquartile range (IQR)
Line in box is usually the median
Whiskers extend to last OBSERVATION within 1 step (usually 1.5*IQR) from end of the box
Any observations beyond whisker are plotted as individual points
value residuals squared residuals
1 12.7 2.88 8.27
2 14.2 4.32 18.7
3 11.3 1.41 1.99
4 7.43 -2.44 5.93
5 8.27 -1.60 2.57
6 12.7 2.81 7.90
7 5.34 -4.53 20.5
8 13.9 4.00 16.0
sum of residuals sum of squares mean
1 -0.0577 359 9.87
For a population the variance \(S_p^2\) is exactly the mean squared distance of the values from the population mean
\[ s{^2}_{pop} = \frac{\sum(x - \bar{x})^2}{N} \]
But this is a biased estimate for the population variance
A biased sample variance will underestimate population variance
n-1 (if you take a large enough sample size, will correct for this)
Small \(s\) = taller, narrower Large \(s\) = squatter, wider
Histograms plot frequency/density of observations within bins
Quantile-Quantile plots plot quantiles of a dataset vs. quantiles of a theoretical (usually normal) distribution
If data is normally distributed, then if we plot the quantiles vs. theoretical normal quantiles (with same n), then we find a linear relationship
The standard deviation tells you how far each score lies from the mean in your datset
Standard error measures the accuracy of a sample mean as an estimate of the population mean.
\[ SE = {s\over \sqrt n} \]
Standard error of the difference measures the accuracy of the difference between two sample means.
\(Mean~difference = \overline x_1 - \overline x_2\)
\(SED = \sqrt{s_1^2\over n_1}+{s_2^2\over n_2}\)
Can be constructed by YOU:
Do you want to know the range within which you will capture the mean…
66% of the time = 1 * S.E.
95% = 1.96 * S.E.
99% = 2.58 * S.E.
The Confidence level (C) is the inverse of the significance level.
\[ \alpha = 1 - C \]
Under realistic scenarios we only ever have an estimate of the standard deviation.
At large enough sample sizes this doesn’t matter, but how large?
The t-distribution is a correction for small sample size and unknown population variance with a normally distributed population
the t-distribution is closely related to a normal distribution, but influenced by sample size n and has fatter tails
Calculate Confidence intervals based on t-distribution
\({95\%~CI} = {\overline x \pm t*SE}\)
Calculating t based on sample size WILL affect confidence intervals
p-values and confidence intervals are both derived from the same underlying theory
a confidence interval that does not include the null hypothesis (e.g. zero difference) has a corresponding p-value
For example a 95% confidence interval that does not include zero is the same as p < 0.05
In this example of a simple linear model, we run the equivalent to a Student’s t-test.
Q. What happens when we have more than two groups in our predictor variable? Why can’t we just do more t-tests?
We have collected data for soil uranium concentrations at three locations on Los Alamos National Lab property: Site A, Site B, and Site C. The data structure is shown here:
A one-way ANOVA can be used to assess whether there is a statistically significant difference in uranium concentration in soil at three locations
What does this one-way/single-factor refer to?
There is a single factor (variable), with at least 3 levels, where we are trying to compare means across the different levels of that factor.
ANOVA does this indirectly by looking at total between and within group variances as a ratio (F).
\[ SSE = \underset{i=1}{n \atop{\sum}}(y_i - \hat{y_i})^2 \]
\[ SSR = \underset{i=1}{n \atop{\sum}}(\hat{y_i} - \overline{y})^2 \]
\[ SST = \underset{i=1}{n \atop{\sum}}(y_i - \overline{y})^2 \]
where:
\(y_i\) = Observed value
\(\hat{y_i}\) = Value estimated by model
\(\overline{y}\) = The Grand Mean
\(SSE + SSR = SST\)
\[ \large F = {SSR / (k-1)\over SSE / (N-k)} = {MSR\over MSE} \]
\(k\) = Total number of groups
\(N\) = numerator degrees of freedom = Total number of observations across all groups
\(N-k\) = denominator degrees of freedom
\(MSR\) = Mean Squares Regression
\(MSE\) = Mean Squares Error
This is a ratio of the between group variance and the the within group variance.
The F-value or ratio of variances, over their respective degrees of freedom will have an F-distribution.
This F distribution is used when we want to compare within and between group variances.
The curve of the distribution depends on the degrees of freedom, and it is always positively skewed
The higher the F-value the greater the signal-to-noise ratio.
For a given value of numerator and denominator degrees of freedom we can look up the probability of observing this ratio under a null hypothesis of identical variances.
If F value is high enough then we might have enough evidence to conclude that samples are likely drawn from populations with different means.
Bill depth - a continuous quantitative variable
Body mass - a continuous quantitative variable
Penguin Species - a categorical variable with three levels
##Dummy variables
In order to represent a categorical variable in a model - we need to convert these cases to ones that are compatible with the maths.
Dummy variables, are coded with 0 and 1’s. With three species we would have the following:
Adelie[0,0] - the reference or intercept species
Chinstrap[1,0]
Gentoo[0,1]
Start by fitting a model that includes your dependent variable(s), covariates and specifies any key interactions
Check the model is a good fit for the observed data
Diagnostic plots of residuals
Multicollinearity
Omitted Variable Bias
Irrelevant variables
Hypothesis testing for model simplification
The next two plots analyze for collinearity and high leverage points.
High Leverage Points are observations that deviate far from the average. These can skew the predictions for linear models.
Collinearity is when features are highly correlated, which can throw off simple regression models. However, when we are using interaction terms - we expect these to be highly correlated with main effects.
If the distributions are skewed, this can indicate problems with the model.
Quantile-Quantile Plot: We can see that several points towards the end of the quantile plot do not fall along the straight-line.
Normal density Plot: The residuals fit a very nice normal distribution
Also known at the Ordinary Least Squares (OLS) model
It aims to fit a regression line which minimises the squared differences between observed and predicted values
The linear regression line is the most likely line given your data if we assume each data point comes from a hypothetical bell curve centered on the regression line
How often do these assumptions really hold true?
A good fit:
How often do these assumptions really hold true?
An ok fit:
How often do these assumptions really hold true?
A poor fit:
Purpose: To check for linearity and homoscedasticity.
Interpretation: Look for a random scatter of points around the horizontal line at y = 0. A funnel-shaped pattern suggests heteroscedasticity.
Purpose: To assess the normality of the residuals
Interpretation: Points should fall along the diagonal line. Deviation from the line indicates non-normality of residuals.
Purpose: To check homoscedasticity and identify outliers
Interpretation: Uses standardised residuals. Constant spread indicates homoscedasticity
Purpose: To identify influential data points(outliers)
Interpretation: Loking for points with high leverage that might affect the regression line. Investigate values > 0.5
The performance package from easystats can produce (among other things) similar plots.
What type of variable is the outcome?
Fit models: check residuals and AICs
janka_ls <- lm(hardness ~ dens, data = janka)
summary(janka_ls)
janka |>
ggplot(aes( x = dens, y = hardness))+
geom_point()+
geom_smooth(method = "lm")
Call:
lm(formula = hardness ~ dens, data = janka)
Residuals:
Min 1Q Median 3Q Max
-338.40 -96.98 -15.71 92.71 625.06
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1160.500 108.580 -10.69 2.07e-12 ***
dens 57.507 2.279 25.24 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 183.1 on 34 degrees of freedom
Multiple R-squared: 0.9493, Adjusted R-squared: 0.9478
F-statistic: 637 on 1 and 34 DF, p-value: < 2.2e-16
The R output for the MASS::boxcox() function plots a maximum likelihood curve (with a 95% confidence interval - drops down as dotted lines) for the best transformation for fitting the dependent variable to the model.
janka_sqrt <- lm(sqrt(hardness) ~ dens, data = janka)
plot(janka_sqrt, which=2)
plot(janka_sqrt, which=3)
ggplot(janka, aes(x = hardness, y = dens)) +
geom_point() + # scatter plot of original data points
geom_smooth(method = "lm", formula = y ~ sqrt(x)) + # regression line
labs(title = "Sqrt Linear Regression with ggplot2",
x = "X", y = "Y") # axis labels and title
janka_log <- lm(log(hardness) ~ dens, data = janka)
plot(janka_log, which=2)
plot(janka_log, which=3)
janka_poly <- lm(log(hardness) ~ poly(dens, 2), data = janka)
plot(janka_poly, which=2)
plot(janka_poly, which=3)
summary(janka_poly)
ggplot(janka, aes(x = hardness, y = dens)) +
geom_point() + # scatter plot of original data points
geom_smooth(method = "lm", formula = (y ~ poly(log(x), 2))) + # regression line
labs(title = "Quadratic Log Linear Regression",
x = "X", y = "Y") # axis labels and title
Call:
lm(formula = log(hardness) ~ poly(dens, 2), data = janka)
Residuals:
Min 1Q Median 3Q Max
-0.22331 -0.05708 -0.01104 0.07500 0.18871
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.1362 0.0168 424.896 < 2e-16 ***
poly(dens, 2)1 3.3862 0.1008 33.603 < 2e-16 ***
poly(dens, 2)2 -0.5395 0.1008 -5.354 6.49e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1008 on 33 degrees of freedom
Multiple R-squared: 0.9723, Adjusted R-squared: 0.9706
F-statistic: 578.9 on 2 and 33 DF, p-value: < 2.2e-16
janka_wls <- lm(sqrt(hardness) ~ dens, weights = 1/sqrt(hardness), data = janka)
plot(janka_wls, which=2)
plot(janka_wls, which=3)
Call:
lm(formula = sqrt(hardness) ~ dens, data = janka, weights = 1/sqrt(hardness))
Weighted Residuals:
Min 1Q Median 3Q Max
-0.55072 -0.19009 -0.03393 0.19082 0.63845
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.01819 0.96647 2.088 0.0443 *
dens 0.76142 0.02201 34.601 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.294 on 34 degrees of freedom
Multiple R-squared: 0.9724, Adjusted R-squared: 0.9716
F-statistic: 1197 on 1 and 34 DF, p-value: < 2.2e-16
janka_glm <- glm(hardness ~ dens, data = janka, family = gaussian(link = "sqrt"))
plot(janka_glm, which=2)
plot(janka_glm, which=3)
A Cullen and Frey plot (skewness-kurtosis)
A graphical view of the distribution of the dependent variable
Reference lines represent expected skew and kurtosis for different distributions
Beta family sits outside of GLM exponential families
\(f(x;\alpha,\beta) = \frac{\beta^\alpha x^{\alpha - 1} e^{-\beta x}}{\Gamma(\alpha)}\)
Shape (α): Determines the shape or curve of the distribution. Higher values of α result in distributions that are more peaked and skewed to the left.
Rate (β): Influences the spread of the distribution. Smaller values of β result in distributions that are more spread out, larger values of β lead to distributions that are less spread out.
Useful for modelling positively skewed continuous non-negative data
Likelihood measures the probability of observing the given data under a specific statistical model, as a function of the model’s parameters.
It represents how well the model explains the observed data.
The log likelihood function in maximum likelihood estimations is usually computationally simpler as it is additive.
A contour plot visualizing the log-likelihood surface for a normal distribution fitted to simulated data.
The contour lines on the plot represent regions of equal log-likelihood values. Darker regions indicate lower log-likelihood values, while lighter regions indicate higher log-likelihood values
The deviance is a key concept in Generalised linear models.
It measures the deviance of the fitted Generalised linear model with respect to a perfect/saturated model for the sample.
Formula:
\(GLM~R^2 = 1 - \frac{D}{D_0}\)
\(linear\_model~R^2= 1 - \frac{SSE}{SST}\)
\(D\) is the deviance of the fitted model.
\(D_0\) is the deviance of the null model (linear model with only an intercept).
The GLM \(R^2\) is not the percentage of variance explained by the model, but rather a ratio indicating how close the fit is to being perfect.
Characteristic |
Beta |
95% CI 1 |
p-value |
|---|---|---|---|
| type | |||
| Control | — | — | |
| Inseminated | 3.4 | -2.6, 9.0 | 0.3 |
| Virgin | -14 | -20, -9.1 | <0.001 |
| thorax | 151 | 130, 171 | <0.001 |
| sleep | 0.02 | -0.08, 0.14 | 0.7 |
| 1
CI = Confidence Interval |
|||
Investigating mayfly abundance in relation to pollution gradient measured in Cumulative Criterion Units (CCU)
Mayflies serve as indicators of stream health due to sensitivity to metal pollution
Higher CCU values indicate increased metal pollution, potentially impacting mayfly presence or absence in the stream
Fitting a linear regression vs a logit-link regression:
The Bernoulli distribution models the probability of success or failure in a single trial of a binary experiment:
where success occurs with probability \(p\)
and failure with probability \({1-p}\)
The odds put our expected values on a 0 to +Inf scale.
The log transformation puts our expected values on a -Inf to +Inf scale.
Now, the expected values can be linearly related to the linear predictor.
For a simple Bernoulli/Binary GLM there are only a few checks that apply:
# Calculate and plot "deviance" residuals
dresid <- resid(mayfly_glm, type = "deviance")
hist(dresid)
# Plot leverage
plot(mayfly_glm, which = 5)
# Plot Cook's distance
plot(mayfly_glm, which = 4)
Residuals in binary glms will never be normally distributed - so a “residuals vs fitted plot” provides little insight. We can try “binned residuals”:
We can use emmeans package to compute estimated marginal means (EMMs) from a fitted GLM, convert them into a tibble format for visualization with ggplot2, overlay the means with uncertainty ribbons and add them on a scatter plot of the original data.
means <- emmeans::emmeans(mayfly_glm,
specs = ~ CCU,
at=list(CCU=c(0:5)),
type='response') |>
as_tibble()
ggplot(Mayflies, aes(x=CCU, y=Occupancy)) + geom_point()+
geom_ribbon(data = means,
aes(x = CCU,
y = prob,
ymin = asymp.LCL,
ymax = asymp.UCL),
alpha = .2)+
geom_line(data = means,
aes(x = CCU,
y = prob)) + # regression line
labs(title = "Logit-link Binomial Regression",
x = "CCU", y = "Occupancy")+
theme_classic(base_size = 14)40:00
Using the malaria dataset fit a binomial glm to look at the predictors of malaria in the Seychelles warbler
Check the model fit
Make predictions
Produce a suitable figure
\(logit(p) = \log \frac{p}{1 - p}\)
\(y_i = Binomial(n,p)\)
It is the collection of Bernoulli trials for the same event
It is represented by the number of Bernoulli trials \(n\)
It is also the probability of an event in each trial \(p\)
This small dataset is from an experiment looking at mortality in batches of flour exposed to different doses of pesticides.
The question of interest is:
How does pesticide dose affect mortality in the flour beetle Tribolium castaneum
# Calculate and plot "deviance" residuals
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
pred <- predict(beetle_glm_weights)
dresid <- resid(beetle_glm_weights, type = "deviance")
# Check for overdispersion
hist(dresid)
# Binned residuals
arm::binnedplot(dresid, pred)
# Plot leverage
plot(beetle_glm_weights, which = 5)
plot(beetle_glm_weights, which = 2)
Standard Error, Confidence intervals and p-values WILL change.
Estimates remain the same
On the afternoon of January 28th, 1986 the space shuttle Challenger experienced a critical failure and broke apart in mid air, resulting in the deaths of all seven crewmembers.
An investigation discovered critical failures in all six of the O-rings in the solid rocket booster.
How did temperature affect the probability of o-ring failure?
What was the probability of all o-rings failing when the temperature was 36 degrees F
40:00
Count or rate data are ubiquitous in the life sciences (e.g number of parasites per microlitre of blood, number of species counted in a particular area). These type of data are discrete and non-negative.
A useful distribution to model abundance data is the Poisson distribution:
a discrete distribution with a single parameter, λ (lambda), which defines both the mean and the variance of the distribution.
In a study by Kilner et al. (1999), the authors studied the begging rate of nestlings in relation to total mass of the brood of reed warbler chicks and cuckoo chicks. The question of interest is:
How does nestling mass affect begging rates between the different species?
This model is inadequate. It is predicting negative begging calls within the range of the observed data, which clearly does not make any sense.
Let us display the model diagnostics plots for this linear model.
Curvature
Funnelling effect
We can look at the kurtosis and skewness of our univariate variable - Begging:
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
fp <- fitdist(cuckoo$Beg, "pois")
fn <- fitdist(cuckoo$Beg, "norm")
fnb<- fitdist(cuckoo$Beg, "nbinom")
plot.legend <- c("normal", "poisson", "negative binomial")
denscomp(list(fn, fp, fnb), legendtext = plot.legend)
qqcomp(list(fn, fp, fnb), legendtext = plot.legend)
Fit a Poisson model to the cuckoo data
Look at the residual plots - how have they changed?
15:00
Call:
glm(formula = Beg ~ Mass + Species + Mass:Species, family = poisson(link = "log"),
data = cuckoo)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.334475 0.227143 1.473 0.14088
Mass 0.094847 0.007261 13.062 < 2e-16 ***
SpeciesWarbler 0.674820 0.259217 2.603 0.00923 **
Mass:SpeciesWarbler -0.021673 0.008389 -2.584 0.00978 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 970.08 on 50 degrees of freedom
Residual deviance: 436.05 on 47 degrees of freedom
AIC: 615.83
Number of Fisher Scoring iterations: 6
VIF, or Variance Inflation Factor, measures the extent of multicollinearity among independent variables in regression analysis.
High multicollinearity, indicated by elevated VIF values, can distort regression coefficients and inflate standard errors, undermining the reliability of the model’s predictions
This dataset contains standardised counts of Monarch butterflies across three different gardens.
Is there a difference in observance of the butterflies between these three sites?
Check the poisson model is a good fit for this data
Get the monarchs.csv dataset
15:00
10:00
Open the parasite dataset and script
Determine whether the poisson model is suitable for this dataset?
The example is an experiment measuring the effect of the parasitic tapeworm Schistocephalus solidus infection on the susceptibility of infection from a second parasite, the trematode Diplostomum pseudospathaceum
The question of interest is:
How does Treatment/Infection status affect trematode counts?
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
fp <- fitdist(parasite$Diplo_intensity, "pois")
fn <- fitdist(parasite$Diplo_intensity, "norm")
fnb<- fitdist(parasite$Diplo_intensity, "nbinom")
plot.legend <- c("normal", "poisson", "negative binomial")
denscomp(list(fn, fp, fnb), legendtext = plot.legend)
qqcomp(list(fn, fp, fnb), legendtext = plot.legend)
library(DHARMa)
plot(simulateResiduals(stick_nb))
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(stick_nb)
| Model | Acronym/Abbr | Source of Overdispersion |
|---|---|---|
| Zero-inflated Poisson | ZIP | Zeros |
| Negative Binomial | NegBin | Positive integers |
| Zero-inflated Negative Binomial | ZINB | Zeros + Positive integers |
20:00
Using the revegetation dataset
Check the model fit
Determine if zero-inflation is required
Validate the model
\({y_i = \beta0_{[i]} + \beta1x+\epsilon}\)
y ~ x + (1|group), data = data)
\({y_i = \beta0_{[i]} + \beta1_{[i]}x+\epsilon}\)
y ~ x + (x|group), data = data)
Checking models with random effects is more complex
check_model will add plots to investigate the distribution of random effects
DHARMa when dealing with increasingly complex models - residual simulations may be preferred
Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
Scaled residual values: 0.036 0.112 0.808 0.22 0.364 0.72 0.656 0.48 0.464 0.172 0.908 0.604 0.228 0.276 0.94 0.428 0.164 0.372 1 0.58 ...